setwd("/mnt/md1200/5/zhaocunyou/wangzhongju/epistasis_analysis_by_wangzhongju/SMR_for_cognitive/data/Davies_MP_2016")
library(readr)

fn =list.files(pattern="txt")
con.file=list.files(pattern = ".csv",path="/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/727case_confile/",full.names=T)
case=read.table("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/727case_confile/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
caseid=gsub("chr","",unique(case$unitID))

col_names=c("con.file","group_file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=8))
names(result)=col_names
result=result[-1,]
for(j in 1:length(fn)){
file=read_delim(fn[j],col_names=T,delim=" ")
file1=file[file$`P-value`<0.05,]
anno=data.frame(Chromosome=unique(file$Chromosome),chr=1:22)
file1=merge(anno,file1,by="Chromosome")
file1$unitID=paste(file1$chr,file1$Position,sep=":")
sigid =as.character(unique(file1$unitID))

for(i in 1:length(con.file)){
con=read.csv(con.file[i],head=T,sep=",")
con$unitID=paste(gsub("chr","",con$Chr),con$Start,sep=":")
conid=unique(con$unitID)
result_tmp=data.frame(matrix(NA,1,ncol=8))
names(result_tmp)=col_names
  result_tmp[,1]=gsub(".hg19_multianno.csv","",(gsub("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/727case_confile//","",con.file[i])))
  result_tmp[,2]=fn[j]
  result_tmp[,3]=length(intersect(caseid,sigid))
  result_tmp[,4]=length(caseid)-length(intersect(caseid,sigid))
  result_tmp[,5]=length(intersect(conid,sigid))
  result_tmp[,6]=length(conid)-length(intersect(conid,sigid))
  result_tmp[,7]=fisher.test(matrix(c(result_tmp[1,3],result_tmp[1,4],result_tmp[1,5],result_tmp[1,6]),nrow = 2))$estimate
  result_tmp[,8]=fisher.test(matrix(c(result_tmp[1,3],result_tmp[1,4],result_tmp[1,5],result_tmp[1,6]),nrow = 2))$p.value
  result=rbind(result,result_tmp)
  print(unique(file$Chromosome))
  print(length(sigid))
}
}


setwd("/mnt/md1200/5/zhaocunyou/wangzhongju/epistasis_analysis_by_wangzhongju/SMR_for_cognitive/data/Davies_NC_2018")
fn =list.files(pattern="txt")
con.file=list.files(pattern = ".csv",path="/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/727case_confile/",full.names=T)
case=read.table("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/727case_confile/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
caseid=gsub("chr","",unique(case$unitID))

col_names=c("con.file","group_file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=8))
names(result)=col_names
result=result[-1,]
for(j in 1:length(fn)){
file=read_delim(fn[j],col_names=T,delim=" ")
file1=file[file$P<0.05,]
file1$unitID=paste(file1$CHR,file1$BP,sep=":")
sigid =as.character(unique(file1$unitID))

for(i in 1:length(con.file)){
con=read.csv(con.file[i],head=T,sep=",")
con$unitID=paste(gsub("chr","",con$Chr),con$Start,sep=":")
conid=unique(con$unitID)
result_tmp=data.frame(matrix(NA,1,ncol=8))
names(result_tmp)=col_names
  result_tmp[,1]=gsub(".hg19_multianno.csv","",(gsub("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/727case_confile//","",con.file[i])))
  result_tmp[,2]=fn[j]
  result_tmp[,3]=length(intersect(caseid,sigid))
  result_tmp[,4]=length(caseid)-length(intersect(caseid,sigid))
  result_tmp[,5]=length(intersect(conid,sigid))
  result_tmp[,6]=length(conid)-length(intersect(conid,sigid))
  result_tmp[,7]=fisher.test(matrix(c(result_tmp[1,3],result_tmp[1,4],result_tmp[1,5],result_tmp[1,6]),nrow = 2))$estimate
  result_tmp[,8]=fisher.test(matrix(c(result_tmp[1,3],result_tmp[1,4],result_tmp[1,5],result_tmp[1,6]),nrow = 2))$p.value
  result=rbind(result,result_tmp)
  print(unique(file$CHR))
  print(length(sigid))
}
}